#include "yocto/math/io/data-set.hpp"
#include "yocto/exception.hpp"
#include "yocto/ios/icstream.hpp"
#include "yocto/sequence/vector.hpp"
#include "yocto/string/vfs-utils.hpp"
#include "yocto/string/conv.hpp"
#include "yocto/ios/ocstream.hpp"
#include "yocto/math/dat/linear.hpp"
#include "yocto/code/utils.hpp"
#include "yocto/math/sig/data-smoother.hpp"

#include <cmath>

using namespace yocto;
using namespace math;

// for overshoot, resampled 1 600

int main( int argc, char *argv[] )
{
    const char *progname = _vfs:: get_base_name(argv[0]);
    try
    {
        if( argc <= 2 )
            throw exception("Usage: %s file.dat dt [tmin filter_window]", progname);
        
        const string filename = argv[1];
        const double dt       = strconv::to_real<double>(argv[2],"dt");
        double       tmin     = 0;
        if( argc>3 )
            tmin = strconv::to_real<double>(argv[3],"tmin");
        
        double filter_window = 0;
        if( argc>4 )
            filter_window = strconv::to_real<double>(argv[4],"filter_window");
        
        vector<double>   X;
        vector<double>   Y;
        data_set<double> ds;
        ds.use(1,X);
        ds.use(2,Y);
        std::cerr << "Loading '" << filename << "'" << std::endl;
        {
            ios::icstream fp( filename );
            ds.load(fp);
        }
        const size_t N = X.size();
        std::cerr << "#Loaded=" << N << std::endl;
        
        //======================================================================
        // put to 0
        //======================================================================
        {
            const double t_org = X[1];
            for( size_t i=1; i <= N; ++i ) X[i] -= t_org;
        }
        const double t_max = X[N];
        std::cerr << "t_max=" << t_max << std::endl;
        
        size_t       Ns = floor(t_max/dt);
        while( (Ns+1)*dt > t_max && Ns>0) --Ns;
        ++Ns;
        std::cerr << "Resampling every   " << dt << std::endl;
        std::cerr << "Producing  #data = " << Ns << std::endl;
        std::cerr << "sampled    max   @ " << Ns*dt << std::endl;
        if( Ns>1 )
        {
            vector<double> Xs(Ns,0);
            vector<double> Ys(Ns,0);
            for( size_t j=1; j <= Ns; ++j )
            {
                const double x = ( Xs[j] = (j-1)*dt );
                Ys[j] = math::linear<double>(x, X, Y);
            }
            
            vector<double> Xf,Yf;
            
           
            for( size_t i=1; i <= Ns; ++i )
            {
                if( Xs[i] >= tmin )
                {
                    Xf.push_back( (Xs[i]-tmin) );
                    Yf.push_back(Ys[i]);
                }
            }
            
            const size_t Nf = Xf.size();
            if( filter_window > 0 )
            {
                std::cerr << "filtering" << std::endl;
                const size_t     np = max_of<size_t>( filter_window/(2*dt), 2 );
                smoother<double> sm(np,np,smoother<double>::uniform,4);
                extender<double> ex( Yf, extension_constant, extension_constant );
                vector<double>   sy( Nf, 0.0 );
                sm.run( sy, ex, 0 );
                for( size_t i=1; i <= Nf; ++i ) Yf[i] = sy[i];
            }
            string outfile = filename;
            _vfs::change_extension(outfile, "resampled.dat" );
            std::cerr << "saving into " << outfile << std::endl;
            {
                ios::ocstream fp( outfile, false );
                for( size_t i=1; i <= Nf; ++i )
                {
                    fp("%.15g %.15g\n", Xf[i]/60.0, Yf[i] ); 
                }
            }
            
        }
        
    }
    catch( const exception &e )
    {
        std::cerr << e.what() << std::endl;
        std::cerr << e.when() << std::endl;
        return -1;
    }
    catch(...)
    {
        std::cerr << "Unhandled exception in " << progname << std::endl;
        return -1;
    }
    
    
    return 0;
    
}
